package cs.ubbcluj.ro.lp;

import java.io.IOException;
import java.util.Iterator;
import java.util.Map;

import cs.ubbcluj.ro.algorithm.AlgorithmBase;
import cs.ubbcluj.ro.algorithm.AlgorithmPredictorCorrector;
import cs.ubbcluj.ro.algorithm.AlgorithmType;
import cs.ubbcluj.ro.algorithm.AlgorithmWeightedCombined;
import cs.ubbcluj.ro.algorithm.AlgorithmPrimalDualSqrt;
import cs.ubbcluj.ro.algorithm.AlgorithmWeighted;
import cs.ubbcluj.ro.algorithm.AlgorithmWeightedSqrt;
import cs.ubbcluj.ro.datatypes.interfaces.Vector;
import cs.ubbcluj.ro.datatypes.interfaces.Matrix;
import cs.ubbcluj.ro.linearsystem.implementations.LinearSystemLU;
import cs.ubbcluj.ro.mps.bridge.implementation.SparseMpsImp;
import cs.ubbcluj.ro.mps.bridge.interfaces.MpsImpInterface;
import cs.ubbcluj.ro.mps.data.ConstrainType;
import cs.ubbcluj.ro.mps.data.interfaces.*;
import cs.ubbcluj.ro.mps.reader.MpsParser;

public strictfp class LP {
	
	Matrix A;
	Vector x;
	Vector y;
	Vector s;
	Vector w;
	Vector b;
	Vector c;
	int rowSize;
	int columnSize;
	int initialSize;
	private Vector l;
	private Vector u;
	private Double optimum_lack;
	
	// constructors
	public LP(String fileName, AlgorithmType algorithmType, double sigma, double rho, double epsilon, double theta) throws IOException {
		
		SparseMpsImp sparseMps = new SparseMpsImp();
		solveProblem(fileName, sparseMps, algorithmType, sigma, rho, epsilon, theta);		
	}
	
	public LP(String fileName, MpsImpInterface mpsInterface, AlgorithmType algorithmType, double sigma, double rho, double epsilon, double theta) throws IOException {
		
		solveProblem(fileName, mpsInterface, algorithmType, sigma, rho, epsilon, theta);		
	}
	
	// the solver functions
	void solveProblem(String fileName, MpsImpInterface mpsInterface, AlgorithmType algorithmType, double sigma, double rho, double epsilon, double theta) {
		try {
			
		    // create the parser
		    final MpsParser parser = new MpsParser(mpsInterface);

		    // parse the file
		    final Mps mps = parser.parse(fileName);		    

		    // number of rows in A
		    this.rowSize = mps.getNumberOfConstrains();

		    // number of columns in A
		    this.columnSize = mps.getCoefficients().getNrOfColumns();

		    // get the A matrix
		    this.A = mps.getCoefficients();

		    // get the B vector
		    this.b = mps.getRightHandSideValues();

		    // get the C vector
		    this.c = mps.getTargetCoefficients();
		    
		    // get the lower bounds vector
		    this.l = mps.getLowerBounds();
		    
		    // get the upper bounds vector
		    this.u = mps.getUpperBounds();
		    
		    this.optimum_lack=0.00;		    
		
		    transformUpperBounds();
		    
			initialPointSettings();
						
			printInitialValues();
			
			AlgorithmBase algorithm = null;
			switch (algorithmType) {
			case PRIMAL_DUAL:
				algorithm = new AlgorithmBase(this, sigma, rho, epsilon, theta);
				break;
			case PRIMAL_DUAL_SQRT:
				algorithm = new AlgorithmPrimalDualSqrt(this, sigma, rho, epsilon, theta);
				break;
			case WEIGHTED:
				algorithm = new AlgorithmWeighted(this, sigma, rho, epsilon, theta);
				break;
			case WEIGHTED_SQRT:
				algorithm = new AlgorithmWeightedSqrt(this, sigma, rho, epsilon, theta);
				break;
			case WEIGHTED_COMBINED:
				algorithm = new AlgorithmWeightedCombined(this, sigma, rho, epsilon, theta);
				break;
			case PREDICTOR_CORRECTOR:
				algorithm = new AlgorithmPredictorCorrector(this, sigma, rho, epsilon, theta);
				break;
			}
			algorithm.solveProblem();
			
		} catch (InterruptedException e) {
			e.printStackTrace();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}
	
	void transformInequalities(Map<String, Constrain> constrains) {
		int iterationIndex = 0;
		Iterator<Map.Entry<String, Constrain>> it = constrains.entrySet().iterator();
	    while (it.hasNext()) {
	        Map.Entry<String, Constrain> constrainEntry = (Map.Entry<String, Constrain>)it.next();
	        Constrain currentConstrain = (Constrain)constrainEntry.getValue();
	        
	        if (currentConstrain.getType() == ConstrainType.L) {	        	
	        	A.addElement(iterationIndex, A.getNrOfColumns(), -1);
	        	c.addElement(c.getNrOfRows(), 0);
	        }
	        
	        if (currentConstrain.getType() == ConstrainType.G) {	        	
	        	A.addElement(iterationIndex, A.getNrOfColumns(), 1);
	        	c.addElement(c.getNrOfRows(), 0);
	        }
	        
	        iterationIndex++;
	    }
	}
	
public void transformLowerBounds(){
		
//		System.out.println("Lower Bounds");
//		l.print();
		int nonZeros=A.getNrOfNonZeros();
		for ( int i=0; i<l.getNrOfRows(); i++ ){
			if ( l.get(i) != 0 ){
				//System.out.println(i);
				for ( int j=0; j<nonZeros; j++){
					if ( A.columnIndex(j) == i ){
						//System.out.println(A.columnIndex(j));
						try {
							
							b.addElement(A.rowIndex(j), b.get(A.rowIndex(j))-A.get(A.rowIndex(j),A.columnIndex(j))*l.get(i));
							
						} catch (InterruptedException e) {
							e.printStackTrace();
						}
					}
				}
				if ( u.get(i) != 0 ){
					try {
						u.set(i, u.get(i)-l.get(i));
					} catch (InterruptedException e) {
						e.printStackTrace();
					}
				}
				optimum_lack +=c.get(i)*l.get(i);
			}
		}
//		System.out.println("b after transforming lower bounds");
//		b.print();
//		
//		System.out.println("optimum lack: "+optimum_lack);
	}
	
	public void transformUpperBounds(){
		
//		System.out.println("Upper Bounds");
//		u.print();	
		int nrOfRows = A.getNrOfRows();
		int nrOfColumns=A.getNrOfColumns();
		int j=-1;
		for ( int i=0; i<u.getNrOfRows(); i++ ){
			if ( u.get(i) != 0 ){
				    j++;
					b.addElement(nrOfRows+j, u.get(i));
					A.addElement(nrOfRows+j, i, 1.0);
					A.addElement(nrOfRows+j, j+nrOfColumns, 1.0);
					c.addElement(nrOfColumns+j,0.0);
			}
		}

//		System.out.println("A after transforming upper bounds");
//		A.print();
//		System.out.println("b after transforming upper bounds");
//		b.print();
//		System.out.println("c after transforming upper bounds");
//		c.print();
	}
	
	public void initialPointSettings() throws InterruptedException {		
		
		// calculate A * A^t (it's needed for further use)
		Matrix AAt = A.product(A.trans());
		
		// calculates the initial y value, given by the formula: A*A^T * y = A * c
		LinearSystemLU linearSystem = new LinearSystemLU();
		Vector y0 = linearSystem.solveSystem(AAt, A.product(c));
		
		// calculates the initial s value, given by the formula: s = c - A^T * y		
		Vector s0 = c.diff(A.trans().product(y0));
		
		// calculates the initial x value, given by the formula: x = A^T * (A * A^T)^-1 * b
		Vector auxX = linearSystem.solveSystem(AAt, b);
		Vector x0 = A.trans().product(auxX);

		// calculate deltaX = max (-1.5 * min(x0), 0)
		double deltaX = x0.min();
		deltaX *= - 1.5;
		if (deltaX < 0) deltaX = 0;
		
		// calculate deltaS = max (-1.5 * min(s0), 0)
		double deltaS = s0.min();
		deltaS *= - 1.5;
		if (deltaS < 0) deltaS = 0;
		
		// auxiliary  variables for the calculation of deltaX0 and deltaS0		
		Vector x0DeltaX = x0.createHelperVector(x0.getNrOfRows());		
		x0DeltaX.set(x0);
		x0DeltaX.sum(deltaX);
		
		Vector s0DeltaS = s0.createHelperVector(s0.getNrOfRows());
		s0DeltaS.set(s0);
		s0DeltaS.sum(deltaS);
		
		// calculate deltaX0 = deltaX + (x0 + deltaX * e)^T * (s0 + deltaS * e) / 2 * sum(s0 + deltaS * e)
		deltaX += x0DeltaX.product(s0DeltaS) / (2 * s0DeltaS.sumOfComponents()); 		
		
		// calculate deltaS0 = deltaS + (x0 + deltaX * e)^T * (s0 + deltaS * e) / 2 * sum(x0 + deltaX * e)
		deltaS += x0DeltaX.product(s0DeltaS) / (2 * x0DeltaX.sumOfComponents());		
		
		// set initial interior point variables
		// y = y0
		y = y0;
		
		// s = s0 + deltaS
		s = s0.sum(deltaS);
		
		// x = x0 + deltaX	
		x = x0.sum(deltaX);	
		
		w = calculateInitialW();
	}
	
	public Vector calculateInitialW() throws InterruptedException {
		Vector xs = x.createHelperVector(x.getNrOfRows());
		xs.set(x);
		xs.productScalar(s);
		xs.sqrt();
		
		return xs;
	}
	
	// getter functions for private variables
	public Matrix getA() {
		return A;
	}
	
	public Vector getX() {
		return x;
	}
	
	public Vector getY() {
		return y;
	}
	
	public Vector getS() {
		return s;
	}
	
	public Vector getW() {
		return w;
	}
	
	public Vector getC() {
		return c;
	}
	
	public Vector getB() {
		return b;
	}
	
	public int getXLength() {
		return x.getNrOfRows();
	}

	/**
	 * rb = A * x - b
	 * @return value of rb
	 * @throws InterruptedException 
	 */
	public Vector getRb() throws InterruptedException {	
		return A.product(x).diff(b);
	}
	
	/**
	 * rc = At * y + s - c
	 * @return value of rc
	 * @throws InterruptedException
	 */
	public Vector getRc() throws InterruptedException {			
		return A.trans().product(y).sum(s).diff(c);		
	}
	
	/**
	 * 
	 * @return matrix given by the formula: A * x * s^-1
	 * @throws InterruptedException
	 */
	public Matrix getAXS_1 () throws InterruptedException {
		Matrix helperMatrix = A.createHelperMatrix(A.getNrOfRows(), A.getNrOfColumns());
		helperMatrix.set(A);
		
		int nrOfColumns = A.getNrOfColumns();
		for (int i = 0; i < nrOfColumns; i++) {
			helperMatrix.scalarProductColumn(i, x.get(i)/s.get(i));			
		}
		
		return helperMatrix;
	}
	
	/**
	 * 
	 * @return matrix given by the formula: A * s^-1
	 * @throws InterruptedException
	 */
	public Matrix getAS_1 () throws InterruptedException {	
		Matrix helperMatrix = A.createHelperMatrix(A.getNrOfRows(), A.getNrOfColumns());
		helperMatrix.set(A);
		
		int nrOfColumns = A.getNrOfColumns();
		for (int i = 0; i < nrOfColumns; i++) {
			helperMatrix.scalarProductColumn(i, 1/s.get(i));
		}
		
		return helperMatrix;
	}
	
	/**
	 * calculates AD2At using the formula: AD2At = A * X *S_1 * At
	 * @return
	 * @throws InterruptedException
	 */
	public Matrix getAD2At() throws InterruptedException {	
		Matrix helperMatrix = A.createHelperMatrix(A.getNrOfRows(), A.getNrOfRows());		
		helperMatrix.set(getAXS_1().product(A.trans()));
		
		return helperMatrix;
	}
	
	
	/**
	 * calculates the value of S^-1
	 * @return
	 * @throws InterruptedException
	 */
	public Vector getInvS() throws InterruptedException {
		return s.inv();
	}
	
	/**
	 * calculates the Primal value
	 */
	public double getPrimal() throws InterruptedException {
		Vector helperVector = x.createHelperVector(x.getNrOfRows());
		helperVector.set(x);
		return helperVector.product(c);
	}
	
	/**
	 * calculates the Dual value
	 */
	public double getDual() throws InterruptedException {	
		Vector helperVector = y.createHelperVector(y.getNrOfRows());
		helperVector.set(y);
		return helperVector.product(b);
	}
	
	// setter functions for private variables
	/**
	 * sets the values of x to the given values
	 */
	public void setX(Vector v) throws InterruptedException {
		this.x.set(v);
	}
	
	/**
	 * sets the values of y to the given values
	 */
	public void setY(Vector v) throws InterruptedException {
		this.y.set(v);
	}
	
	/**
	 * sets the values of s to the given values
	 */
	public void setS(Vector v) throws InterruptedException {
		this.s.set(v);
	}
	
	/**
	 * sets the value of w to the given values
	 * @param v
	 * @throws InterruptedException
	 */
	public void setW(Vector v) throws InterruptedException {
		this.w.set(v);
	}
	
	void printInitialValues() {
//		System.out.println("Initial values");
//		System.out.println("A = ");
//		A.print();
//
//		System.out.println("\nb = ");
//		b.print();
//
//		System.out.println("\nc = ");
//		c.print();

//		System.out.println("\nx = ");
//		x.print();

//		System.out.println("\ny = ");
//		y.print();

//		System.out.println("\ns = ");
//		s.print();
		
//		System.out.println("\nw = ");
//		w.print();
	}
	
}
